SUBROUTINE set_near_lattice
USE param
USE global
IMPLICIT NONE

INTERFACE
   FUNCTION reallocint(p,n)
     USE param
     IMPLICIT NONE
     INTEGER,POINTER,DIMENSION(:) :: reallocint,p
     INTEGER,INTENT(in) :: n
   END FUNCTION reallocint
END INTERFACE

INTEGER :: i,j,k,l,m,n,p,q
INTEGER :: cx,cy,cz,cell_no,neighbor_cell_no,lattice_no,belong_to
REAL(DBL) :: dist,dist2,px,py,pz,dx,dy,dz
INTEGER :: n_vac,n_sia,n_nor,n_odd,n_defect,n_bound,atom_no1,atom_no2,d1,d2,d3,total_atom

INTEGER,POINTER,DIMENSION(:) :: defect_list,defect_type_list
INTEGER :: n_alloc_defect_list,alloc_size

1001 FORMAT(3X,'>>> Setting nearest atoms for lattice...')
1002 FORMAT('Filename: ',A17,1X,'V: ',I6,3X,'I: ',I6,3X,'Odd: ',I6,3X,'D: ',I6)
1003 FORMAT(I6,1X,'V1',1X,1PE13.6,1X,1PE13.6,1X,1PE13.6)
1004 FORMAT(I6,1X,'I1',1X,1PE13.6,1X,1PE13.6,1X,1PE13.6,1X,I2,1X,I2,1X,I2)
1005 FORMAT(3X,'PKA_E',1X,1PE13.6,1X,'eV',1X,I6)
2001 FORMAT('US FE')
2002 FORMAT(3X,I6)
2003 FORMAT('C',1X,1PE13.6,1X,1PE13.6,1X,1PE13.6)
2004 FORMAT('I',1X,1PE13.6,1X,1PE13.6,1X,1PE13.6)
2005 FORMAT('Fe',1X,1PE13.6,1X,1PE13.6,1X,1PE13.6)

PRINT 1001

alloc_size=100
ALLOCATE(defect_list(alloc_size))
ALLOCATE(defect_type_list(alloc_size))
n_alloc_defect_list=alloc_size

! --- Check number of atoms
dummy_line1=9
DO i=1,n_lammps_file
   ! Initialize lattice
   DO j=1,n_lattice
      lattice(j)%n_near_atom=0
   END DO
   ! Read from source file and Count point defects
   OPEN(UNIT=1,FILE=lammps_file(i)%source,STATUS='OLD')
   ! Read dummy line
   DO j=1,dummy_line1
      READ(1,'(A)') eatline
   END DO
   ! Read atoms
   DO j=1,n_atom
      READ(1,'(A)') eatline
      READ(eatline,*,END=101) (eatword(k),k=1,mn_word)
101   CONTINUE
      ! label no
      read(eatword(1),*) p
      ! Read x position
      READ(eatword(3),*) px
      ! Read y position
      READ(eatword(4),*) py
      ! Read z position
      READ(eatword(5),*) pz
      atom(j)%no=p
      atom(j)%px=px
      atom(j)%py=py
      atom(j)%pz=pz
      ! Find the cell box to which belong
      cx=FLOOR(px/dist_a0_ang)
      cy=FLOOR(py/dist_a0_ang)
      cz=FLOOR(pz/dist_a0_ang)
      cell_no=cell_box_map(cx,cy,cz)
      if (cx.lt.0) px=px+box_width
      if (cx.ge.nside_int) px=px-box_width
      if (cy.lt.0) py=py+box_width
      if (cy.ge.nside_int) py=py-box_width
      if (cz.lt.0) pz=pz+box_width
      if (cz.ge.nside_int) pz=pz-box_width
      dist=1.0E+08
      belong_to=0
      DO k=1,cell_box(cell_no)%n_neighbor
         neighbor_cell_no=cell_box(cell_no)%neighbor_cell_list(k)
         DO l=1,cell_box(neighbor_cell_no)%n_lattice
            lattice_no=cell_box(neighbor_cell_no)%lattice_list(l)
            dist2=(lattice(lattice_no)%px+box_width*cell_box(cell_no)%bx(k)-px)**2+&
                  (lattice(lattice_no)%py+box_width*cell_box(cell_no)%by(k)-py)**2+&
                  (lattice(lattice_no)%pz+box_width*cell_box(cell_no)%bz(k)-pz)**2
!            dist2=(lattice(lattice_no)%px-px)**2+&
!                  (lattice(lattice_no)%py-py)**2+&
!                  (lattice(lattice_no)%pz-pz)**2
!            print *,dist2
            IF(dist.GT.dist2) THEN
               dist=dist2
               belong_to=lattice_no
            END IF
         END DO
      END DO
!      print *,belong_to,dist
!      pause
      IF(belong_to.EQ.0) STOP 'ERR:Null lattice number.'
      ! Count the number of atoms belong to the site
      lattice(belong_to)%n_near_atom=lattice(belong_to)%n_near_atom+1
      lattice(belong_to)%near_atom_list=>&
           reallocint(lattice(belong_to)%near_atom_list,lattice(belong_to)%n_near_atom)
      lattice(belong_to)%near_atom_list(lattice(belong_to)%n_near_atom)=j
   END DO
   CLOSE(1)

   ! Prepare to write
   ! Count point defect and save defect list
   n_bound=0
   n_vac=0
   n_sia=0
   n_nor=0
   n_odd=0
   n_defect=0
   total_atom=0
   DO j=1,n_lattice
      IF(lattice(j)%removed) THEN
         n_bound=n_bound+1
         CYCLE
      END IF
      IF(lattice(j)%n_near_atom.EQ.0) THEN
         ! Vacancy
         n_vac=n_vac+1
         n_defect=n_defect+1
         total_atom=total_atom+1
         IF(n_alloc_defect_list.LT.n_defect) THEN
            n_alloc_defect_list=n_alloc_defect_list+alloc_size
            defect_list=>reallocint(defect_list,n_alloc_defect_list)
            defect_type_list=>reallocint(defect_type_list,n_alloc_defect_list)
         END IF
         defect_list(n_defect)=j
         defect_type_list(n_defect)=0
      ELSE IF(lattice(j)%n_near_atom.EQ.1) THEN
         n_nor=n_nor+1
      ELSE IF(lattice(j)%n_near_atom.EQ.2) THEN
         ! SIA
         n_sia=n_sia+1
         n_defect=n_defect+1
         total_atom=total_atom+2
         IF(n_alloc_defect_list.LT.n_defect) THEN
            n_alloc_defect_list=n_alloc_defect_list+alloc_size
            defect_list=>reallocint(defect_list,n_alloc_defect_list)
            defect_type_list=>reallocint(defect_type_list,n_alloc_defect_list)
         END IF
         defect_list(n_defect)=j
         defect_type_list(n_defect)=2
      ELSE
!         cx=FLOOR(lattice(j)%px/dist_a0_ang)
!         cy=FLOOR(lattice(j)%py/dist_a0_ang)
!         cz=FLOOR(lattice(j)%pz/dist_a0_ang)
!         cell_no=cell_box_map(cx,cy,cz)
!         print '("Lattice info:",i0,1x,i0,1x,1pe10.3,1x,1pe10.3,1x,1pe10.3)'&
!              ,j,cell_no,lattice(j)%px,lattice(j)%py,lattice(j)%pz
!         do p=1,lattice(j)%n_near_atom
!            q=lattice(j)%near_atom_list(p)
!            cx=FLOOR(atom(q)%px/dist_a0_ang)
!            cy=FLOOR(atom(q)%py/dist_a0_ang)
!            cz=FLOOR(atom(q)%pz/dist_a0_ang)
!            cell_no=cell_box_map(cx,cy,cz)
!            print '("Atom:",i0,1x,i0,1x,i0,1x,1pe10.3,1x,1pe10.3,1x,1pe10.3)'&
!                 ,q,atom(q)%no,cell_no,atom(q)%px,atom(q)%py,atom(q)%pz
!         end do
!         pause
         n_odd=n_odd+1
         n_defect=n_defect+1
         total_atom=total_atom+lattice(j)%n_near_atom
         IF(n_alloc_defect_list.LT.n_defect) THEN
            n_alloc_defect_list=n_alloc_defect_list+alloc_size
            defect_list=>reallocint(defect_list,n_alloc_defect_list)
            defect_type_list=>reallocint(defect_type_list,n_alloc_defect_list)
         END IF
         defect_list(n_defect)=j
         defect_type_list(n_defect)=3
      END IF
   END DO

   ! Test for counting
   PRINT 1002,TRIM(ADJUSTL(lammps_file(i)%source)),n_vac,n_sia,n_odd,n_defect
!   IF(n_vac.NE.n_sia) PRINT *,'WARNING: n_vac is not equal to n_sia.'
!   IF(n_atom.NE.n_bound+n_vac+n_nor+n_sia+n_odd) STOP 'ERR: Counting is invalid.'
!   IF(n_odd.GT.0) PRINT *,'WARNING: n_odd is larger than 0.'

   ! Write to out file
   OPEN(UNIT=1,FILE=lammps_file(i)%outfile,STATUS='REPLACE')
   WRITE(1,1005) pkae,n_defect
   DO j=1,n_defect
      IF(defect_type_list(j).EQ.0) THEN
         ! Vacancy
         WRITE(1,1003) j,&
                       lattice(defect_list(j))%px*aton,&
                       lattice(defect_list(j))%py*aton,&
                       lattice(defect_list(j))%pz*aton
      ELSE
         ! SIA or odd type
         ! Setting the direction of align of SIA
         atom_no1=lattice(defect_list(j))%near_atom_list(1)
         atom_no2=lattice(defect_list(j))%near_atom_list(2)
         dx=atom(atom_no1)%px-atom(atom_no2)%px
         dy=atom(atom_no1)%py-atom(atom_no2)%py
         dz=atom(atom_no1)%pz-atom(atom_no2)%pz
         IF(dx.GT.0) THEN
            d1=1
         ELSE IF(dx.LT.0) THEN
            d1=-1
         ELSE
            d1=0
         END IF
         IF(dy.GT.0) THEN
            d2=1
         ELSE IF(dy.LT.0) THEN
            d2=-1
         ELSE
            d2=0
         END IF
         IF(dz.GT.0) THEN
            d3=1
         ELSE IF(dz.LT.0) THEN
            d3=-1
         ELSE
            d3=0
         END IF
         WRITE(1,1004) j,&
                       lattice(defect_list(j))%px*aton,&
                       lattice(defect_list(j))%py*aton,&
                       lattice(defect_list(j))%pz*aton,&
                       d1,d2,d3
      END IF
   END DO
   CLOSE(1)

   ! Write to xmakemol file
   OPEN(UNIT=1,FILE=xmfn,STATUS='OLD',ACCESS='APPEND')
   WRITE(1,2002) total_atom+n_xml_corner
   WRITE(1,2001)
   DO j=1,n_xml_corner
      WRITE(1,2003) xml_corner(j)%px,xml_corner(j)%py,xml_corner(j)%pz
   END DO
   DO j=1,n_defect
      IF(defect_type_list(j).EQ.2) THEN
         atom_no1=lattice(defect_list(j))%near_atom_list(1)
         atom_no2=lattice(defect_list(j))%near_atom_list(2)
         WRITE(1,2004) atom(atom_no1)%px,atom(atom_no1)%py,atom(atom_no1)%pz
         WRITE(1,2004) atom(atom_no2)%px,atom(atom_no2)%py,atom(atom_no2)%pz
      END IF
   END DO
   DO j=1,n_defect
      IF(defect_type_list(j).EQ.0) THEN
         WRITE(1,2005) lattice(defect_list(j))%px,&
                       lattice(defect_list(j))%py,&
                       lattice(defect_list(j))%pz
      END IF
   END DO
   DO j=1,n_defect
      IF(defect_type_list(j).EQ.3) THEN
         DO k=1,lattice(defect_list(j))%n_near_atom
            atom_no1=lattice(defect_list(j))%near_atom_list(k)
            WRITE(1,2004) atom(atom_no1)%px,atom(atom_no1)%py,atom(atom_no1)%pz
         END DO
      END IF
   END DO
   CLOSE(1)

END DO

DEALLOCATE(defect_list)
DEALLOCATE(defect_type_list)

END SUBROUTINE set_near_lattice
